Introduction

Patterned distributions of individual plant species and species groups are likely to reflect patterns in environmental conditions and constraints imposed on plant survival, growth, and reproductive ability by abiotic and biotic factors, including physiological tolerance of levels of various environmental conditions (e.g., soil moisture, [Ca], wind speed), competition with other plants, and interactions with herbivores and pathogens. If variation in the abundance of focal species characterized by different growth forms or physiologies is correlated with environmental conditions, and not due merely to limited plant dispersal, it can suggest that those forms or physiologies are adapted to different states or levels of those conditions.

Sarracenia is a genus of 11 species of carnivorous plants, occurring on the coastal plains of the eastern United States, that produce pitcher-shaped leaves filled with liquid that function as pitfall traps.

Givnish (personal communication) observed a patchy distribution of S. alata in the extensive wet, sandy pine savannas at Buttercup Flats near Wiggins, Mississippi in the DeSoto National Forest. Such wet pine savannas of the Gulf Coastal plain of the United States contain an extraordinary diversity and abundance of carnivorous plants occurring down slope from mesic longleaf pine (Pinus palustris Mill.) dominated uplands (Brewer, 1999). On these sites, S. alata tends to co-occur with crayfish mounds made of clay material (Givnish, pers. obs.), suggesting that the plant’s locally patchy distribution may have originated from a hidden variability in the physical composition of the substrate (Brewer, 1999). Specifically, the presence of clay lenses in the subsoil may create localized anoxic zones that prevent nitrogen fixation by soil microbes, leading to wet, N-poor microsites that favor the establishment of carnivorous plants, due to low levels of soil nitrate and hostile, anoxic conditions that can work against root uptake of soil nutrients generally. I therefore aim to test the hypothesis that scattered but hidden clay lenses help generate the patchy distribution of Sarracenia, creating zones that are poor in nitrate and dissolved oxygen.

pattern1

pattern1

Thirteen transects crossing the described vegetation pattern were set up randomly within in the savanna at Buttercup Flats facing a compass bearing of 160 degrees. Each transect contained a pitcher-absent microsite in the middle and 10-m of pitcher-present microsites on each end. Transects 11, 12, and 13 intersected multiple pitcher-absent patches and were therefore evaluated at more locations compared to the other transects. Soil and tissue nutrients, ground Coverage of all species, and soil clay content clay content were measured at these sites using a stratified random design.

Analysis of Soil Nutrients

Plant Root Simulators (PRS, available from Western Ag) with cation and anion exchange resins be planted in pairs of three in two locations along the transect: 1) a random location at the pitcher-present microsite in the beginning ten meters of the transect and 2) a random location in each pitcher-absent microsite contained within the transect. Probe burial time was 46-47 days. All soil nutrients were measured in μg/cm2.

Each nutrient is be evaluated independently for its association with microsite type (pitcher plant-present, pitcher plant-absent).

soilNPK = read.csv(file = "soilNPK.csv", header = TRUE)
soilNPK$Transect = factor(soilNPK$Transect) # interpreted as factor, not numerical values
str(soilNPK)
'data.frame':   31 obs. of  6 variables:
 $ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ NO3      : num  0.6 5.2 0.9 0.9 3.8 2.8 3.2 1.5 1 1.7 ...
 $ NH4      : num  2.3 18.9 1.6 2.6 3.1 1.9 7.3 2 2.8 2.3 ...
 $ K        : num  128 321 105 130 147 ...
 $ P        : num  0.9 0.3 0.2 0.3 0.3 0.3 0.4 0.5 0.3 0.9 ...

Soil Potassium

data visualization

layout(matrix(1:2,1,2))
hist(soilNPK$K, main="", col="tan", breaks=20)
qqnorm(soilNPK$K, main = '', ylab="soil K")
qqline(soilNPK$K)

The histogran for soil potassium concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. Three outliers (large K values) can be seen in the Q-Q plots. they are on 3 different transects, but all in the “-” microsite: pitcher plants absent (note that one of these transects, #13, has another observation for “-” microsite)

subset(soilNPK, K>200 )
   Microsite Transect NO3  NH4      K   P
2          -        1 5.2 18.9 321.01 0.3
8          -        4 1.5  2.0 229.97 0.5
31         -       13 2.2  3.8 225.55 0.6

Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transect 1 is in red, transect 4 is in blue).

transect.color = rep("black",13)
transect.color[c(1,4)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, K, legend=FALSE, col=transect.color))
plot(K~Microsite, data=soilNPK)

The outliers are much “less” outliers when we consider the log(potassium) values:

layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(K), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(K), main = '', ylab="soil log(K)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(K), legend=FALSE, col=transect.color))
plot(log(K)~Microsite, data=soilNPK)

From this, 4 transects have K higher in the pitcher-present area, and 9 transects have K higher in the pitcher-absent area.

statistical analysis

If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.

At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.

A new vector is created specifically for the paired t-test.

table(soilNPK[1:2])
         Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
        - 1 1 1 1 1 1 1 1 1  1  4  2  2
        + 1 1 1 1 1 1 1 1 1  1  1  1  1
logKmeans = with(soilNPK, tapply(log(K), Microsite:Transect, mean))
K_minus = logKmeans[1:13]
K_plus = logKmeans[14:26]
K_minus  # log values, in fact
     -:1      -:2      -:3      -:4      -:5      -:6      -:7      -:8 
5.771472 4.868534 4.652435 5.437949 4.959482 3.810876 4.453417 4.704925 
     -:9     -:10     -:11     -:12     -:13 
4.279994 4.211387 4.529404 4.798584 4.838847 
K_plus
     +:1      +:2      +:3      +:4      +:5      +:6      +:7      +:8 
4.850858 4.654817 4.988935 4.259718 4.596129 4.636960 4.755915 4.425086 
     +:9     +:10     +:11     +:12     +:13 
4.076859 4.863913 4.371218 4.717427 4.620847 

A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.

t.test(K_plus, K_minus, paired = TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -0.74329, df = 12, p-value = 0.4716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4531952  0.2226375
sample estimates:
mean of the differences 
             -0.1152788 

K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.47).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.23).

Kmeans = with(soilNPK, tapply(K, Microsite:Transect, mean))
K_minus = Kmeans[1:13]
K_plus = Kmeans[14:26]
t.test(K_plus, K_minus, paired=TRUE)

    Paired t-test

data:  K_plus and K_minus
t = -1.2574, df = 12, p-value = 0.2325
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -72.50872  19.44179
sample estimates:
mean of the differences 
              -26.53346 

A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.

The linear model is constructed with K (or log(K)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between K and microsite (p=0.22).

fit <- lm(K ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: K
          Df Sum Sq Mean Sq F value Pr(>F)
Transect  12  39671  3306.0  1.0959 0.4207
Microsite  1   4940  4939.6  1.6375 0.2179
Residuals 17  51282  3016.6               
drop1(fit, test="F")
Single term deletions

Model:
K ~ Transect + Microsite
          Df Sum of Sq   RSS    AIC F value Pr(>F)
<none>                 51282 257.75               
Transect  12     40938 92220 251.94  1.1309 0.3981
Microsite  1      4940 56222 258.60  1.6375 0.2179
layout(matrix(1:4,1,4, byrow=T))
plot(fit)

We can re-do this analysis, using log(K) values this time. Still no evidence of an association (p=0.47). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of K (untransformed).

fit <- lm(log(K) ~ Transect + Microsite, data=soilNPK)
anova(fit)
Analysis of Variance Table

Response: log(K)
          Df  Sum Sq Mean Sq F value Pr(>F)
Transect  12 2.12012 0.17668  0.9749 0.5064
Microsite  1 0.10048 0.10048  0.5544 0.4667
Residuals 17 3.08095 0.18123               
drop1(fit, test="F")
Single term deletions

Model:
log(K) ~ Transect + Microsite
          Df Sum of Sq    RSS     AIC F value Pr(>F)
<none>                 3.0810 -43.571               
Transect  12   2.15191 5.2329 -51.150  0.9895 0.4954
Microsite  1   0.10048 3.1814 -44.576  0.5544 0.4667
plot(fit)

Soil Nitrogren

Soil nitrogen is evaluated as both Total Nitrogen and its components:
Soil Total Nitrogen is the sum of soil nitrate (NO3) and soil ammonium (NH4).

soilNPK$TN <- soilNPK$NO3 + soilNPK$NH4

Analysis of Sedge tissue nutrients

Aboveground biomass of 6 species of sedges in the Rhynchospora genus (R. megalocarpa, R. compressa, R. gracilenta, R. plumosa, R. rariflora, R. harveyi var harveyi) were collected, dried, and analyzed for % N, P, and K content, as a measure of the availability of these nutrients to non-carnivorous plants. All species are C3.

All tissue nutrients were measured in % of Aboveground Dry Mass.

tissuenutrients <- read.table("tissueNPK.csv", header = TRUE, sep=",")
tissuenutrients$Transect = factor(tissuenutrients$Transect)
str(tissuenutrients)
'data.frame':   29 obs. of  5 variables:
 $ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
 $ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ TN       : num  1.55 1.07 0.63 0.7 1.09 0.81 1.01 0.52 1.08 0.55 ...
 $ P        : num  0.14 0.0869 0.08 0.06 0.0868 ...
 $ K        : num  0.67 0.59 0.81 0.61 0.52 0.89 0.6 0.63 0.61 0.9 ...

NPK Limitation Diagram

N:P, N:K, and P:K ratios from tissue nutrient measurements are compared to critical ratios put forth by Olde Venterink et al. (2003) to determine whether plants are P or P + N limited (N:P > 14.5, K:P >3.4), K or K + N limited (N:K > 3.1, K:P < 3.4), or N limited (N:P < 14.5, N:K < 2.1). The relationship between N:P:K ratios and type of nutrient limitation is illustrated in the subsequent triaxial graph.

library(ggtern)
plot <- ggtern(data = tissuenutrients, aes(x = TN, y = P, z = K)) +
  geom_Tisoprop(value=.3226, color='darkgreen') + # T=Top apex: K/(K+N)=.3226 i.e N/K=1/.3226-1=2.1
  geom_Tisoprop(value=.2439, color='lightgreen') + # Top apex: K/(K+N)=.2439 i.e N/K=1/.2439-1=3.1
  geom_Risoprop(value=.06452, color='darkblue') + # Right apex P/(P+N)=0.145: N/P=1/.06452-1=14.5
  geom_Lisoprop(value=.7727273, color='darkred') + # Left apex K/(K+P)=: P/K=1/0.7727273-1=, K/P=1/0.2941=3.4
  geom_point(aes(fill = Microsite),
             size = 2,
             shape = 21,
             color = "black") +
  ggtitle("NPK limitation diagram") +
  labs(fill = "Microsite") +
  theme_rgbw() +
  theme(legend.position = c(0,1),
        legend.justification = c(1, 1))
plot

Critical ratios are represented by vertical lines.

Plants in both pitcher-present and pitcher-absent microsites are most strongly N-limited (?). CA: I would say that they are most strongly P-limited.

About N limitation:

  • All sites have N:P higher than 14.5 (blue isoproportion line). (so none are N limited, correct?).
  • None of the N:K values are higher than 3.1 (light green line), but some are higher than 2.1 (darkgreen line). (these would definitely not be N limited, correct?)

About P limitation: all points are super far from the top apex (100% P), close to the horizontal line across that apex, where there is 0% P. More specifically:

  • all sites have K:P above 3.4 (red line), and
  • all sites have N:P above 4.5 (blue line).

The points are overlapping a lot around the line N:K = 2.1, so it’s hard to see if there are mostly pitcher-present or pitcher-absent sites above/below this line. But let’s tabulate: we get that there are only 4 sites with TN:K above 2.1. Of these 4 sites, 3 are in a pitcher-absent patch (2 along the same transect 11, the third in transect 13), and the 4th is in a pitcher-present area. It goes in the expected direction, but 3 versus 1 is not statistically different from a half/half ratio.

with(tissuenutrients, table(TN/K<2.1, Microsite))
       Microsite
         -  +
  FALSE  3  1
  TRUE  13 12
subset(tissuenutrients, TN/K>2.1)
   Microsite Transect   TN        P    K
1          +        1 1.55 0.140000 0.67
24         -       11 1.12 0.075815 0.50
25         -       11 1.18 0.075889 0.55
29         -       13 1.17 0.098511 0.48
transect.color = rep("black",13)
transect.color[c(2,13,11)] = c("red", "blue","purple")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, TN/K, legend=F, col=transect.color))
plot(TN/K ~ Microsite, data=tissuenutrients)

From this we see that in 10 of the 13 transect, TN:K is higher in the pitcher-present area than in the pitcher-absent patch. Does that run counter to expectations? (in the plot above and those below, transect 2 is in red, transect 13 is in blue, and transect 11 is in purple.)

If Phosphorous is limiting, we could analyze the ratios P/xxx. For example:

layout(matrix(1:4,2,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/K, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/TN, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(K+TN), legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(P+K+TN), legend=F, col=transect.color))

We if look at the proportion of P, i.e. P/(P+K+TN), we still don’t find any evidence of association with microsite (pitcher plant present/absent).

fit = lm(P/(P+K+TN) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")
Single term deletions

Model:
P/(P + K + TN) ~ Transect + Microsite
          Df  Sum of Sq        RSS     AIC F value  Pr(>F)  
<none>                  0.00095217 -271.40                  
Transect  12 0.00156373 0.00251590 -267.22  2.0528 0.09461 .
Microsite  1 0.00019328 0.00114545 -268.04  3.0448 0.10144  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
layout(matrix(1:4,1,4))
plot(fit)

How about N+P limitation?

layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/K, legend=T, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/(TN+P+K), legend=F, col=transect.color))

fit = lm((TN+P)/(TN+P+K) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")
Single term deletions

Model:
(TN + P)/(TN + P + K) ~ Transect + Microsite
          Df Sum of Sq      RSS     AIC F value Pr(>F)
<none>                 0.086167 -140.74               
Transect  12  0.087773 0.173940 -144.37  1.2733 0.3246
Microsite  1  0.005863 0.092030 -140.84  1.0206 0.3284
layout(matrix(1:4,1,4))
plot(fit)

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X Yosemite 10.10.5

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggtern_2.2.0  ggplot2_2.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.11        knitr_1.15.1        magrittr_1.5       
 [4] MASS_7.3-45         munsell_0.4.3       colorspace_1.3-2   
 [7] lattice_0.20-34     stringr_1.1.0       plyr_1.8.4         
[10] tools_3.3.2         grid_3.3.2          gtable_0.2.0       
[13] htmltools_0.3.5     yaml_2.1.14         lazyeval_0.2.0     
[16] rprojroot_1.1       digest_0.6.11       assertthat_0.1     
[19] bayesm_3.0-2        tibble_1.2          latex2exp_0.4.0    
[22] energy_1.7-0        tensorA_0.36        gridExtra_2.2.1    
[25] compositions_1.40-1 robustbase_0.92-7   evaluate_0.10      
[28] rmarkdown_1.3       labeling_0.3        stringi_1.1.2      
[31] DEoptimR_1.0-8      scales_0.4.1        backports_1.0.4    
[34] boot_1.3-18         proto_1.0.0